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ABSTRACT 

A time-dependent model for pair creation in a pulsar magnetosphere is devel¬ 
oped. It is argued that the parallel electric field that develops in a charge-starved 
region (a gap) of a pulsar magnetosphere oscillates with large amplitude. Elec¬ 
trons and positrons are accelerated periodically and the amplitude of the oscilla¬ 
tions is assumed large enough to cause creation of upgoing and downgoing pairs at 
different phases of the oscillation. With a charge-starved initial condition, we find 
that the oscillations result in bursts of pair creation in which the pair density rises 
exponentially with time. The pair density saturates at N± ~ Eg/(87rm e c 2 r thr ), 
where Eq is the parallel electric field in the charge-starved initial state, and r t hr 
is the Lorentz factor for effective pair creation. The frequency of oscillations 
following the pair creation burst is given roughly by u osc = eE 0 /(8m e cT t ^ r ). A 
positive feedback keeps the system stable, such that the average pair creation 
rate balances the loss rate due to pairs escaping the magnetosphere. 


1. Introduction 

A key issue in pulsar and black hole electrodynamics is the injection of plasma on open 
magnetic field lines. In the absence of adequate plasma, referred to as a ‘gap’ where the 
system is said to be ‘charge-starved’, the rotation of the compact object induces a component 
of the electric field parallel to magnetic field lines, E\\, which accelerates any test particle to 
very high energy. Such ‘primary’ particles reach energies at which the photons they emit, 
through curvature radiation (CR) or inverse Compton scattering (ICS), exceed the threshold 
for effective pair creation, leading to a pair cascade that populates the magnetosphere with 
‘secondary’ pairs. The charge-starved E\\ is screened once the pair plasma is sufficiently 
dense to provide the Goldreich-Julian (GJ) charge density. To maintain the electric charge 
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density close to the GJ value, the plasma must be injected continuously at an average rate 
which equals the loss rate due to escape of charged particles along open magnetic lines. 

In a conventional class of pulsar models, (e.g., the reviews by Michel 1991, Beskin, 
Gurevich and Istomin 1993, Mestcl 1998), which we refer to as gap-plus-PFF models, the 
system is assumed to be in a steady state in the corotating frame. The gap, the region of pair 
creation and the region where E\\ is screened are assumed to be spatially separated. Charge 
starving in an ‘inner gap’ near the stellar surface results in pair creation that is strongly 
concentrated in a narrow range of heights, referred to as a pair formation front (PFF) or 
a pair production front. Reflecting secondaries (positrons if the primaries are electrons) 
provide a net, quasi-static charge density such that the PFF may be regarded as a thin 
surfaces with a surface charge density. The resulting additional electrostatic field screens 
out the initial E\\ above the PFF (e.g., Fawley et al. 1977; Arons and Scharlemann 1979; 
Shibata et al. 1998). In some pulsar models pairs are also assumed to be produced in an 
‘outer gap’ (Cheng et al. 1986), and current conservation can be achieved by the currents 
through the inner and outer gaps being part of a global circuit (Shibata 1991). More recent 
models differ quantitatively, rather than qualitatively, from earlier models in that the Lorentz 
factors of the secondary pairs are smaller (Zhang and Harding 2000, Hibschman and Arons 
2000, Arendt and Eilek 2002) and the pair creation is over a more extensive region, with 
the heights of the inner and outer gaps perhaps overlapping (Shibata et al. 2002). Despite 
extensive studies of such models, it remains unclear how and where the radio emission is 
generated. One observational estimate of the location of the source of the radio emission 
(Blaskiewicz,Cordes and Wasserman 1991; Gangadhara and Gupta 2001) suggests that it 
is between the putative inner and outer gaps. Although the possible high-energy emission 
processes are well understood, there is no consensus on the location of the high energy 
emission region, with strong arguments for both the inner gap (Harding and Muslimov 
1998) and the outer gap (Romani 1986). The lack of success of pulsar models in explaining 
the observed emissions (e.g., Beskin 1999) leaves open the possibility that one or more of the 
basic assumptions in the model may be incorrect. A critical review of the assumptions was 
given by Michel (2004), who favored a radically different model involving a large Coulomb 
held associated with charge-separated regions of the magnetosphere (Krause-Polstorff and 
Michel 1985). 

In this paper we relax the steady-state assumption, arguing that stability of models 
involving steady-state pair creation at a PFF is questionable. In fact early models that 
did not make the steady-state assumption pre-date later models that did. In particular, 
Sturrock (1971) argued that a steady state is not possible and that the pairs escape in a 
sequence of sheets, and Ruderman and Sutherland (1975) invoked pair creation in bursts, 
called sparking. One specific difficulty with a steady-state model relates to the reflection 



of secondaries at the PFF. Suppose E« accelerates electrons upward, so that it accelerates 
secondary positrons downward. In a steady-state model most of the secondary positrons 
are assumed to continue propagating outward, requiring that the screening of E\\ above 
the PFF be sufficient to prevent their reflection. However, when this condition is imposed 
explicitly (Shibata et al. 1988), it requires such implausibly large densities that it seems 
untenable. If a substantial fraction of the positrons are reflected, such positrons become 
downward propagating counterparts of the primary electrons, and by the same sequence of 
processes can, under some conditions, create a second PFF near the stellar surface (Harding 
and Muslimov 1998). Reflection of the created electrons at this second PFF would in turn 
greatly increase the number of primary electrons, leading to a rapidly increasing rate of pair 
creation. Such a runaway pair creation would not just short out the initial E y, but overshoot 
it, setting up oscillations with a large amplitude, of order the initial E\\. Fine tuning is 
required in order to have sufficient number of reflected positrons to allow screening, but 
not too many to cause an overshoot setting up oscillations. The stability of the PFF to 
oscillations cannot be explored within the framework of conventional gap-plus-PFF models 
because of the steady-state assumption: the screening is postulated to be independent of 
time. 

Here we explore an alternative model in which E\\ is oscillatory, with an amplitude of or¬ 
der the unscreened E \\. There is some correspondence between spatially localized phenomena 
in steady-state gap-plus-PFF models and temporally localized phenomena in an oscillatory 
model. The inner gap is replaced by an initially charge-starved region that we also refer 
to as a gap; pair creation associated with the oscillations leads to a burst of pair creation 
sufficient to screen the initial E\\. The spatially localized PFF in a gap-plus-PFF model 
is replaced by temporally localized pair creation near the phases where \E\ \ is maximum. 
The important new feature in an oscillatory model, is that the electric field is determined 
primarily by inductive effects, as opposed to electrostatic effects in a gap-plus-PFF model. 
In fact, we assume that the Goldreich-Julian charge density is small compared with the 
instantaneous density of electrons or positrons, such that it can be neglected to a first ap¬ 
proximation and included through a small asymmetry in the oscillatory model. There is a 
steady-state current (along a given magnetic flux tube) that is a free parameter in either 
type of model. However, by hypothesis there is no time-varying current in a steady-state 
model, whereas the time-varying current is intimately related to the oscillating electric field 
through the induction equation in the model discussed here. 

The model is introduced in section 2, and numerical results for some illustrative cases 
are presented in section 3. The model is discussed further in section 4. 
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2. Oscillating pair-creating electric field 


In a frame corotating with the star Maxwell’s equations can be written in the form, 


V-E = 

4vr(p-p G j), 

(i) 

V B = 

o, 

(2) 


10B 

(3) 

V x E = 

c dt ’ 


4tt 1 dE 

(4) 

V x B = 

(j J.r) r\, 1 

c c at 


where pgj is the GJ charge density, and the term is a combination of the fields and their 
derivatives, and is given explicitly in (Fawley et al. 1977). The electromagnetic fields, current 
density, and charge density are in the corotating frame. As shown below, if the difference 
electric charge density initially satisfies \p— pcj| ~ IpgjIj then the amplitude of the oscillating 
current j is of the order of cp G j or larger. In this case we find j a/j ~ QR/c -C 1, and can 
therefore neglect the term j/j in eq. (4). 

Here we construct a simple ID, time-dependent model, and consider purely parallel 
oscillations only. We define s to be the distance along a magnetic field line, and denote 
by E\\ and Ej_ the field-aligned (parallel) and cross-field (perpendicular) components of the 
electric field, respectively. We ignore drift motion (both inertial and ExB) of charged 
particles, and assume the motion of pairs to be restricted to the direction along magnetic 
field lines. We further neglect variations of the electric field and charge density across field 
lines; that is, we assume that Eu and p are functions of s and t only. Equations (l)-(4) then 
reduce to 

4vr(p e -p G j), (5) 

-4vr(j|| - jo), (6) 

0. (7) 


dE \\ 

ds 

dE\\ 

~dt 

V x E = 


Here j 0 = b • (V x B) is time independent, where b is a unit vector along the magnetic field. 
Our neglect of variations across magnetic field lines implies in addition djo/ds = 0. The 
constant jo is a free parameter in our model, associated with the global DC current flowing 
along the corresponding magnetic field line. Equations (5) and (6) are related by virtue of 
charge conservation. Using V • j + dp e /dt = 0, and the fact that j 0 is divergence free, we 
obtain: 


d_ 

dt 



0 . 


( 8 ) 
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Consequently, Poisson’s equation in the form (5) is automatically satisfied, and is needed 
only to determine the initial electric held E\\(t = 0, s) for a given choice of initial charge 
distribution. 

In steady-state models j\\ = j 0 is implicit at all times, so that dE\\/dt = 0. Equation 
(5) is then solved numerically together with the equation of motion for the accelerating 
pairs, and some prescription for pair creation is adopted. In most cases it has been shown 
that for reasonable pair multiplicities (number of secondary pairs per primary particles) the 
model predicts some returning positrons, which are usually ignored. However, any change in 
the number of returning positrons, as necessarily occurs in any perturbation away from the 
steady state, implies a change in the current density inside the gap from the assumed value 
jo- According to eq. (6), any nonzero difference j\\ — jo leads to a temporally varying E\\. 
For a steady-sate model to be stable, the resulting varying E\\ must eventually be damped 
so that the steady-state is restored. However, pair creation should lead to an increase in 
|j|| — jo |, until the system overshoots, leading to oscillations in E\\. Moreover, it is not 
at all clear that for arbitrary initial conditions, the pair cascade process can relax to the 
putative steady-state. Below we demonstrate that indeed a plausible mode of operation is 
large amplitude oscillations that control the pair creation rate via a positive feedback, with 
no indication that the oscillations damp, and hence no indication that the system converges 
to a steady state. 

We adopt a treatment in which the plasma is modeled as a two component fluid, con¬ 
sisting of electrons and positrons. We End it convenient to introduce the equations in a 
covariant formalism (Greek indices run over 0,1,2,3, with signature of the metric +2) and 
then to revert to the one-dimensional (along the parallel axis) form for our specific calcu¬ 
lations. The processes incorporated in the model are emission of curvature radiation (CR) 
photons and one-photon creation of pairs on the external magnetic field. 

Let 

Tjf = h±n±U±U± + p±g al3 (9) 

denote the stress energy tensor of electrons (—) and positrons (+) with n±, U±, p± and h± 
their proper density, 4-velocity, partial pressure and specific enthalpy, respectively, and let 
Q/2 denote the total creation rate per unit volume of electrons (or positrons). The continuity 
equations are 

d a {n±Ul) = Q/2. (10) 

The energy and momentum equations can be expressed as, 

3«C Tf) = ±en ± F^ a Ul - + Q'l (11) 

where F af3 = -F 0a is the Maxwell 4-tensor. The terms S ± and Q /3 account, respectively, for 
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radiative losses due to the emission of CR photons, and the change in energy and momentum 
of the pairs due to pair creation. 

The 4-current density is given by 

f = e(n + U» -n-UH). (12) 

The charge conservation relation = 0 is implied by eq. (10). 

The projection of eq. (11) on the 4-velocity U >3 yields an equation for the change of 
specific entropy, cr±, of each fluid: 

n±T±U±d a cr± = Siu ±0 - Qlu± g - h±Q/ 2, (13) 

where T± are the corresponding temperatures. The first term on the RHS of eq. (13) accounts 
for the decrease in specific entropy due to radiative losses, and the last two terms account 
for the change in the specific entropy resulting from conversion of photons to pairs. 

Using eqs (9), (10), (11), and (13), we obtain the equation of motion of each fluid: 

n±h±U±d a U± = ±en±F\Ul + {-SI + Q%)Gl a - Qfd aP± , (14) 

where 0((' 3 = g al3 + U±U±. 

In order to solve the set of equations we need to make a number of simplifying assump¬ 
tions, and in doing so we seek the simplest possible model that retains the features that we 
consider essential. 

We make the following simplifications. First, under the conditions considered below we 
anticipate the pressure forces to be much smaller than the electric force. We therefore neglect 
the last term on the RHS of eq. (14). Second, we suppose that the specific entropy of each 
fluid is roughly conserved along streamlines, that is, we take U±d a a± = 0 in eq. (13). We 
emphasize that intense pair creation does not necessarily imply that the entropy per particle 
must increase. In fact, under the conditions anticipated, the average energy of CR photons is 
much smaller than the energy of radiating pairs, and so we don’t expect the specific entropy 
to change much during most of the pair creation burst, except for some short initial phase. 
At any rate, the details should not affect dramatically the general conclusions drawn below. 
Third, we ignore any ExB drift, and consider only acceleration along magnetic held lines. 
Fourth, we adopt simple models for the source terms and that are intended to model 
the effects of CR; we ignore ICS which is mode difficult to model. 

We adopt a local coordinate frame in which F 1 0 = E\\. We denote by U = U 1 the 
component of 4-velocity along the magnetic held, by v the corresponding 3-velocity, by 
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N± = n±T± the pair density as measured in the laboratory frame, and by = S±/N±, 
q]_ = Q±/N±, and q± = Q/2N±, respectively, the radiative drag per particle, momentum 
gain per particle due to photon conversion, and pair multiplicity rate. Under the above 
assumptions equations (14) and (10) reduce to: 


dU± 

~dt 


dU+ 


v±- 


ds 


dN± d^ T , 

~ r + its { n±v±) 


± £V 

h± 

Qf 2- 


h± h± 


(15) 

(16) 


The enthalpies h± may be treated as free parameters of the model (they are assumed to be 
time independent ). The maximum Lorentz factor, j max , that a test particle can acquire 
during the initial pair creation burst is limited essentially by radiative drag. To estimate 
7max during the pair creation burst, we equate the electric force with the CR loss per unit 
length to obtain: 7 max ~ 2 x 10 5 £ , ^' /i . For a pulsar having B < B c: where B c = m 2 c’/eh 2 is 
the critical field strength, the vacuum electric held near the surface is E s ~ Q S B S R S < 5 x 10 8 
V/cm, yielding 7 max < 10' . This implies that the characteristic energy of CR photons is, 
e 7 ± = 1.57j), ftyhc/P b < 10 5 m e c 2 , where ps ~ R s is the radius of curvature of the magnetic 
held lines along which the pairs are accelerated. Hence, the energy of newly created pairs 
is much smaller than the maximum energy m e c 2 7 max . Consequently, the fraction of pairs 
having maximum Lorentz factor is expected to be rather small during the pair creation 
burst. Indeed, as shown below, during the pair creation burst the average Lorentz factor is 
limited by mass loading, and is typically well below 7 max - The amplitude of oscillations of 
the 4-velocity decreases as the pair density increases, until it reaches a value at which the 
energy of CR photons is just below the threshold for a single photon annihilation, at which 
point pair creation ceases. 


A detailed derivation of the energy distribution of pairs involves solving the relevant 
kinetic equations, and a simpler approach is needed for the present purpose. We assume 
that the spread in Lorentz factors of particles in each fluid is small with respect to the bulk 
Lorentz factor, and take the energy loss term to be = 2e 2 cT\_/?>p 2 B for CR, where T_ (T + ) 
is the bulk Lorentz factor of the electron (positron) fluid. Linder the reasonable assumption 
that the particles in each fluid are distributed isotropically in the fluid rest frame, we obtain 



s±v± 


2e 2 cT±U± 
3 Pi 


(V) 


The characteristic energy of a CR photon is e 7 ± = 1.5T 8 he/ps — 10” lo r T 8 (ps/10 6 cm) eV. 
Provided the energy of the emitted photon satisfies e 7 sin -0 > 2 m e c 2 , where is the angle 
between the photon momentum and the direction of the local magnetic field, the photon can 



decay into an electron-positron pair. The lifetime for this decay is assumed to be short: all 
photons are assumed to decay instantaneously once the energy of the particle emitting the 
CR photons exceeds an effective threshold value, T t h r . The number of photons of character¬ 
istic energy produced per unit time per unit volume by the electron fluid is approximately 
d 2 Ny/dtdV ~ N_s°_ /e 7 _ ~ 10 2 ALT_, and likewise for the positron fluid. These photons 
are converted into pairs over a distance As <C p B if T± S> T t hr- Thus, the source terms 
associated with pair creation can be approximated as 

Q = a + N + + a-N_, (18) 

q]_ = (a + iV + e 7+ + a_iV_e 7 _)c/2A-i = , (19) 

where a± = 10 2 T± for T± > T t h r , and a± = 0 otherwise. The threshold value above which 
CR photons annihilate inside the gap is a function of position, but for simplicity we take it 
to be a fixed value (except for one example, as explained below), which is a free parameter 
in our model. The value we adopt in our calculations is T t hr = 10 6 . 

Once pairs are accelerated to Lorentz factors in excess of T t h r , the pair density rises 
rapidly. The bulk Lorentz factor grows until the right hand side of eq. (15) vanishes, at 
which point it saturates. Using eqs (17) and (19) we find that this happens roughly when 
eE\\c 2 = h±q±U±. When the pair system becomes roughly symmetric, the Lorentz factors of 
the fluids saturate at a value given by 

r + = T_ ~ lO 4 ^ 72 , (20) 

where Eu is defined below. Our interpretation is that the bulk Lorentz factors are limited by 
mass loading. Specifically, as the bulk Lorentz factors of the accelerating fluids approach this 
value, the pair creation rate increase to a level at which the electric potential energy goes to 
create additional pairs that move at the same bulk Lorentz factor rather than accelerating 
the fluids. 

The above equations can be rendered dimensionless by employing the normalization: 
E = eLE/m e c 2 , N = N/N B , j 0 = j 0 /(eN B c), and measuring velocities in units of c, distances 
in units of L and time in units of u^ 1 . Here N B = Q s B s /2irec, u B = Aite 2 N B /im e) and 
L = c/lv b , with B s being the magnetic field strength at the stellar surface, and = 2 nP is 
the angular velocity of the star. Note that in these units the vacuum electric field near the 
surface is given as 

E s = eLD, s B s R s /m e c 3 = 3 x lO e (B s /B c ) 1/2 (P/l s _1 ) 1/2 (i2 s /10 6 cm). (21) 

In what follows we use this normalization with B s = B c , R s = 10 6 cm, and P = 1 s -1 . 
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3. Results 

To obtain our numerical results we integrated eqs (6), (12), (15), (16), (17), (18), and 
(19) on a uniform computational grid with 10 4 cells. In all the examples presented below 
the inner boundary is at the stellar surface, and the outer boundary above the pair creation 
region. The boundary condition imposed at the outer boundary of the computational domain 
is the zero-gradient conditions, which minimizes the effect of the boundary condition on the 
solution. We made some runs with other boundary conditions and found little differences in 
most cases. 

In the first example that we considered we made the simplest assumptions to illustrate 
the essential features of the model. In this model all quantities are assumed to be initially 
homogeneous within the gap, and the pair creation rate is also homogeneous with a threshold 
Tthr = 10 6 . This is intended to model sparking in a region close to a charge starved stellar 
surface. It is readily seen that under these conditions the system remains homogeneous at all 
times, so that the variations are purely temporal. An example is shown in fig. 1, where the 
electric held (normalized to the initial value E 0 = E\\(t = 0)), and 4-velocity, normalized to 
the threshold Lorentz factor (T t hr = 10 6 in this example) are plotted against dimensionless 
time, for two different initial conditions. The upper two panels correspond to initial electric 
held E 0 = —10 4 , and the lower panels to E 0 = —10 6 . The evolution of the corresponding 
electron density is exhibited in hg. 2; the positron density equals the electron density. The 
initial pair densities in both cases were taken to be N± = 10 -2 , and the initial velocities 
U- = 0.95 and U + = 0. The parameter j 0 was set equal to the initial current density, viz., 
jo = j\\(t = 0) = —6.9 x 10 ~ 3 N b c. We assumed cold fluids: h± = m e c 2 . 

As clearly seen from hg. 1, the 4-velocity increase rapidly initially until it reaches the 
value given by eq. (20), at which point a burst of pair creation commences, as can be 
seen from hg. 2. This is followed by a phase during which the pair density and, hence, the 
electric current density, j'm , rises exponentially while the bulk velocity remains approximately 
constant. As ju grows, the rate of change of the electric held, given by eq. (6), increases. 
The electric held eventually passes through zero and reverses sign leading to a deceleration 
of the bulk how. The system then continues to oscillate. As seen from hg. 1, the amplitude 
of oscillations of the 4-velocity damps over several periods until it reaches the pair creation 
threshold, T thr , at which point the pair cascade process terminates abruptly. 

The sawtooth shape of the electric held oscillations is a consequence of a relativistic 
effect. The electric current density is proportional to the 3-velocity of the charged particles, 
which is approximately constant at ±c during each cycle, except for a short period during 
which the particles become non-relativistic briefly, and the current density changes sign 
abruptly. To see this more clearly, consider the oscillations after the pair creation burst 
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ceases and the pair density saturates. Then N + ~ AT_ = No is independent of time, and so 
j| = e(N + v + — N_v_) ~ 2 eN 0 V- is also time independent as long as U± 3> 1. Eq. (6) then 
implies that Eu changes linearly with time during most of the oscillation period. Specifically, 


f i/r- for 

11 ~ 0 \ (T_ + T+ - t)/T+ for T_ < f < T_ + 2T + , 


( 22 ) 


where T± = u)^Eq/(N ± jo), with E 0 being the amplitude of the electric field oscillations, 
roughly equals the strength of the initial field, and N = N + = AL. The sign change of the 

_1 /q 

electric current occurs over a fraction T thr / of the period. When j 0 ^ 0 the oscillations are 
asymmetric in time. In the example shown in fig. 1, where jo is negative, the time it takes 
the electric field to change from Eu = —E 0 to E\\ = +E 0 is slightly longer than the time it 
takes it to change back. As a consequence the time averaged 4-velocities, ([/_) and (U + ), 
do not vanish as in the case of symmetric oscillations. In fact eqs (6), (12) and (15) imply 
ec(N + (U + ) — N_(U_)) = jo as expected. This asymmetry is small in this example, where the 
DC current jo is much smaller than the oscillating current jn, and is not well resolved in fig. 1. 
By employing eqs (6) and (12) and ignoring this small asymmetry we estimate the oscillation 
period of the saturated system (after pair creation effectively ceases) to be T = u^iAEo/N). 
Integrating eq. (15) from U — 0 to U — U max , using the above result for the period T, and 
taking the amplitude of the 4-velocity oscillations to be U max ~ T thr , yields the saturated 
pair density in terms of the initial electric field strength and the threshold Lorentz factor: 


N ~ N P 


E n 


El 


2 T thr 87T m e c 2 T thl ' 


(23) 


Substituting the last equation into the above expression for the oscillation period we obtain 


T = LU 


-isr 


thr 


B 


8 m e c 

—^1 

eEo 


thr- 


(24) 


We find these analytic estimates to be in excellent agreement with the numerical results. 
Note that there is a positive feedback on the system; any loss of plasma from the gap, as 
expected in more realistic situations in which the gap is finite, would result in a reduction 
of the pair density N and a consequent increase of U max . This would lead to regeneration 
of plasma inside the gap to compensate for the losses. The average pair creation rate would 
then be equal to the loss rate. 


As a second example we take the initial electric field to be the vacuum held of an aligned 
rotator near the axis. The vacuum held decreases with height above the stellar surface, and 
we take this spatial variation into account by adopting a simple model for it. To be more 
specihc, we take E\\(s,t = 0) = E s ( 1 + s/R s )~ 4 , with E s given by eq. (21). Because the 
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magnetic field declines with height (increasing r), the absorption length of a single photon 
at a given energy increases with height, in a manner that depends on the geometry of the 
magnetic field. For illustration we take the source terms Q and q 1 as in eqs (18) and (19), 
but with a± = 10 2 r ± exp[—(s/i? s ) 2 ] if r± > r thr and a± = 0 otherwise. In the cases shown 
in figs 3-5, the initial pair density is N+(s, t = 0) = N-(s, t — 0) = 10 -2 for s < R s and zero 
for s > R s , such that initial charge density is p e (s,t = 0) = e(N + — iV_) = 0. The initial 
4-velocities were chosen such that the initial electric current density is equal to j 0 . We found 
that the final results are insensitive to the choice of the initial densities and velocities of 
electrons and positrons, provided they are much smaller than the final values. The evolution 
of the electric field is shown in fig. 3, where the electric field is plotted as a function of distance 
from the surface (located at s = 0). As seen from this time sequence, a wave pattern is formed 
above the surface, and propagates outwards. This is a consequence of the inhomogeneity of 
the system; since, as discussed above, the frequency of saturated parallel oscillations T~ l is 
roughly proportional to the local electric field strength E\\(s,t = 0), regions closer to the 
stellar surface, where the electric field is stronger, oscillate faster, as clearly seen in fig. 3. 
The corresponding 4-velocity is displayed in fig. 4: the amplitude approaches the (local) 
threshold value T t h r , as in the first example. The amplitude of the electric field oscillations 
near the surface become smaller with time. At c Ost = 600 the electric charge density reduces 
by a factor of about 10 2 relative to its value at c ost = 400, although the pair density does 
not change much, cf. fig. 5. As seen from fig. 4, at c ost = 600 the oscillations of 4-velocity 
near the stellar surface become asymmetric, and it appears that both type of charges are 
convected towards the surface, although we emphasize that at this stage a better resolution 
is needed to resolve the structure near the surface. 

Unfortunately, our numerical resolution is insufficient to follow the system for times 
longer than shown, and it is not clear what the system evolves towards after very long 
time. There is an indication in our calculations that after a sufficiently long time the system 
becomes chaotic. Physically, we expect that after pair creation effectively ceases, the rate 
of pair creation falls below the replacement rate needed to overcome escape of pairs, the 
oscillations break up and a charged-starved configuration is re-established. The evolution 
described above then repeats. However, we cannot confirm this expectation within the 
framework of the present model, which does not include the escape of particles explicitly. 

In our final example, we take the initial gap structure to be as in the model considered 
by Shibata et al. (1998). To be more precise, the initial electric field is given by Eu(s,t = 
0) = s(l — s/s 0 ) for s < s 0 and E\\ = 0 for s > s 0 , cf- fig. 6. In the example shown in 
fig. 6, So = 0.05i? s . For this choice of parameters the initial electric field is much smaller 
than the vacuum field given by eq. (21), but still large enough to accelerate pairs well above 
r thr . Since so <C R s the variation of the magnetic field is expected to be rather small and 
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we, therefore, take pair creation to be uniform in the gap. We find similar results as in 
the previous examples. Oscillations are set up locally inside the gap, and a wave pattern 
develops with time by virtue of the non-uniformity of the initial conditions. Again we are 
unable to follow the system for times longer than shown. 


4. Discussion and conclusions 

The preliminary calculations described in the preceding section demonstrate that pair 
creation in pulsar and black hole magnetospheres is likely to be strongly oscillatory in nature 
rather than quasi-stationary. The essential difference between the oscillatory model proposed 
here and stationary models is in the nature of the response to the E\\ that unavoidably 
develops in a pulsar magnetosphere. In an oscillatory model, this response is essentially 
inductive, intimately coupled to an oscillating current through the induction equation. In 
a stationary model En is essentially electrostatic and the induction equation is irrelevant 
by hypothesis. Our analysis shows that a partially unscreened E\\ tends to develop large 
amplitude oscillations whenever the electric current deviates locally from the global DC 
current flowing in the system. Such deviations seem unavoidable in realistic situations. This 
suggests to us that quasi-stationary models are unstable to temporal perturbations, and that 
the natural state of pair production in a pulsar magnetosphere is oscillatory. Although in 
this paper we focus on pulsars, we expect similar behavior in magnetospheres of Kerr black 
holes. 

Our results show that given an initial charge-starved configuration, the oscillations start 
with a brief episode during which the pair density rises exponentially, resulting from copious 
pair creation at phases of the oscillation around the maximum values of Ay |. The resulting 
charge (and mass) loading of the oscillating system gives rise to an increase in the frequency of 
oscillations (at a fixed height), and a decrease in the amplitude of the 4-velocity oscillations. 
The pair density saturates when the maximum 4-velocity, which is initially well above the 
threshold required for pair creation, approaches this threshold. The subsequent oscillations 
are maintained at a level such that the net creation rate of pairs in the gap balances escape. 
However, the limitations of our analysis precludes explicit treatment of the escape in the 
simplest form of our model. The oscillations in the saturated state are relatively stable, due 
to a positive feedback whereby any reduction in the pair density leads to an increase in the 
Lorentz factor of oscillating pairs and, hence, to a net increase in the pair creation rate. The 
saturated pair density and corresponding frequency of oscillations depend only on the initial 
value of the parallel electric field, and the threshold Lorentz factor for pair creation. 

In situations where the initial state is not spatially uniform, the frequency of oscillations 



varies with position in the gap. Our calculations show a wave-like pattern developing in the 
inhomogeneous case, and it appears that the wave-like oscillations propagate to regions 
outside the gap. Our limited resolution does not allow us to follow the evolution of the 
system for long enough time to examine how the oscillations are transmitted to distant 
regions. 

We regard the calculations presented in this paper as providing a basis for an alterna¬ 
tive oscillatory model to the gap-plus-PFF model for the electrodynamics of pulsar magne¬ 
tospheres. We make a number of major simplifications that need to be relaxed in a more 
realistic model, but we expect the main features identified here to remain in more realistic 
models. We comment specifically on three of our simplifying assumptions. First, the global 
DC current is a free parameter in our model, and the average charge density is determined 
from initial conditions. This means that the global structure of the magnetosphere is not 
addressed within the framework of our analysis, only the microphysics of charge-starved 
regions. Second, we model the pair creation in a highly idealistic ways: pairs are created 
only by particles above a specified threshold, r t h r , and the pair creation is assumed to be 
instantaneous. This crude approximation should be adequate to describe the effects of an 
initial burst of pair creation, starting from a charge-starved configuration, but more accurate 
approximations are needed to describe the effects of ongoing pair creation once the amplitude 
of the oscillations approaches a steady state. Also, the simplified model of pair creation is 
derived assuming that the photons are due to CR, and we do not claim that it is realistic 
for resonant ICS (Sturner 1995, Luo 1996), which is thought to dominate in some pulsars 
(Harding and Muslimov 1998). Third, our assumption that there is no spread in Lorentz 
factors (the electrons and positrons are cold in their respective instantaneous rest frames) 
may not be justified under conditions such that the spread in Lorentz factors associated with 
the newly created pairs is comparable to the maximum Lorentz factor due to the oscillating 
electric held. This is particularly relevant to in the early phase when a burst of pair creation 
occurred. Furthermore, our neglect of the initial velocity of newly created pairs during the 
initial burst is justified only in case of sufficiently strong electric fields, for which the accel¬ 
eration length is much smaller than the characteristic size of the system. The inclusion of 
such effects requires the use of kinetic theory, and this greatly complicates the analysis. 

The possible observational implications of the model are not as straightforward as might 
at first appear. Superficially, one might expect the large-amplitude oscillations to be reflected 
in temporal structures in the observed emission. For example, one could argue that the period 
of the oscillations is plausibly of the same order as the time scale of the variations observed in 
the micro-structure in the radio bursts. However, the emission pattern is strongly dependent 
on relativistic beaming effects, and the underlying oscillation frequency cannot be readily 
deconvolved from the Doppler effect and relativistic beaming. A specific prediction of the 
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model arises from the approximate symmetry between upward and downward propagating 
particles: this suggests that one might observe emission from downward propagating particles 
above the pole on the opposite side of the compact object to the observer. However, in order 
to see either high-energy from the opposite hemisphere, its ray path must not intersect the 
compact object. To see radio emission, not only must it miss the star, but also it must be 
able to propagate through the region of closed held lines around the star where it can be 
absorbed or scattered through a variety of processes. The flux of particles impinging on the 
star must heat the polar cap region, and this can lead to observable thermal X-ray emission. 
In this context the model is similar to the model of Harding and Muslimov (1998) with two 
PFFs. A careful analysis of the predictions of the model are needed before one can test it 
with more conventional gap-plus-PFF models, for which predictions are also uncertain. 

The oscillatory model raises some interesting possibilities in connection with the radio 
emission mechanism. Two mechanisms favored by the model are linear acceleration emission 
(Melrose 1978; Rowe 1995) associated with the large-amplitude oscillations, and plasma-type 
emission at the phase where the relative flow of the electrons and positrons is nonrelativistic 
and Langmuir waves can grow and be converted into escaping radiation (Weatherall 1994). 
These and other possible radio emission mechanisms require a more detailed discussion than 
is appropriate here. 

AL acknowledges discussions with V. Usov, and support by an ISF grant for the Israeli 
Center for High Energy Astrophysics. 
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Fig. 1.— Evolution of the E\\ (left panels) and the 4-velocity (right panels) in a large- 
amplitude oscillation with a uniform initial electric held E$ = —10 4 (upper panels) and 
E 0 = —10 6 (lower panels). The 4-velocity saturates at U/T th = 1, below which pair creation 
becomes ineffective, and the system oscillates after the charge screening overshoots. 
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Fig. 2.— Evolution of the pair density, showing an exponential increase until the 
overshoots, after which it remains roughly constant as the system oscillates. The d. 
corresDonds to the case with En = —10 4 in fie. 1. and the solid line to En = —10 6 . 
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Fig. 3.— Evolution of the system starting from a charge-starved condition with the vacuum 
electric held decreasing as illustrated in the top left panel. As time increases, oscillations 
are set up with their frequency and amplitude decreasing with height, as shown. 


















